Time-convolutionless stochastic Schrodinger equation for open quantum systems 
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Stochastic methods form a general framework for investigating the generally non-Markovian dy- 
namics of a quantum-mechanical system coupled to an environment. They promise to be com- 
putationally superior to the master-equation approach, which is numerically expensive for large 
dimensions of the system Hilbert space. A stochastic Schrodinger equation that is local in time but 
nevertheless reproduces the dynamics of a non-Markovian master equation would be particularly 
useful. Here, we derive such an equation, which can be solved with moderate numerical cost, com- 
parable to that of a Markovian system. We also discuss a portable algorithm for the generation of 
the noise associated with the non-Markovian dynamics. 



The use of stochastic methods to investigate the dy- 
namics of a physical system coupled to an external bath 
has a long history dating back to Einstein and Langevin 
[l|, [2[ . The idea behind this approach is that the many 
degrees of freedom of the bath induce random motion 
in the system [3|47[. Classically, this is due to collisions 
between the particles of the bath and of the system and 
can be described by a Langevin equation for certain vari- 
ables of the system. Quantum- mechanically, the random- 
ness is introduced by transitions between different states 
of the system induced by the environment and can be 
described by a stochastic Schrodinger equation (SSE). 
Alternatively, one can derive statistical descriptions av- 
eraging over many realization of the stochastic process, 
leading to the Fokker-Planck equation for the distribu- 
tion function of a classical system and the master equa- 
tion for the reduced density operator of a quantum sys- 
tem |3H7| . Assuming the equivalence between the master 
equation and the SSE, the latter has sometimes been seen 
as a "quick and dirty" way to obtain the solution of the 
former, especially when accurate results were not neces- 
sary. Recently, however, it has been shown that the two 
approaches are not equivalent and the master-equation 
description might fail when the Hamiltonian is stochas- 
tic or depends on the instantaneous state of the system 
|8( . The numerical solution of the master equation also 
scales poorly with the number of states kept in the cal- 
culation since it is an equation of motion for a matrix in 
state space, whereas a Schrodinger equation is an equa- 
tion for a vector. 

The dynamics of an open system is generally non- 
Markovian. However, the solution of a non-Markovian 
master equation is difficult because it involves the eval- 
uation of a convolution integral which depends on the 
history of the system. Therefore, one often employs a 
Markov approximation, which replaces the memory ker- 
nel in this convolution integral by a S function. We are 
interested in the thermal relaxation dynamics: In contact 
with an external bath, the system will approach an equi- 



librium state. It can be shown that in the non-Markovian 
case there exists an exact condition that the bath corre- 
lation function has to satisfy for the system to reach ther- 
mal equilibrium [9( . This condition is no longer satisfied 
if the Markov approximation is made [9( . 

In order to study non-Markovian dynamics it would 
be advantageous to have a SSE that is local in time but 
is nevertheless able to reproduce the dynamics induced 
by a non-Markovian master equation. In this Letter, we 
introduce a local-in-time (time-convolutionlcss) SSE and 
show its equivalence to a non-Markovian master equa- 
tion. Finally, we illustrate this SSE by studying the 
thermal relaxation of an electronic system coupled to the 
electromagnetic field in a cavity. 

Our starting point is a standard second-order master 
equation [3|, |a|. This non-Markovian master equation 
(NMME) is quite generic, no restrictions on the form 
of the coupling between the system and the environ- 
ment are imposed, except that this coupling is bilinear, 
-Hint = ^^2 a S a <8> B a , in the operators S a and B a from 
the system and the bath, respectively. Here we assume 
that 5*0, and B a are hermitian operators, the extension 
to the more general case where only if; nt is hermitian is 
straightforward. Under the assumptions of weak system- 
bath interaction and factorization of the full density op- 
erator at the initial time t = 0, the equation of motion 
for the reduced density operator p of the system is given 
by (setting ft = 1) 



dp(t) 
dt 



\[H, p(t)] + A 2 £&> ^aC*) - M a (t)} (1) 



with 



M a (t) = J2 / dTc ab (t,T)e-^^-^S bP (T)e^^-^ 



(2) 
In this NMME, H is the Hamiltonian of the sys- 
tem and the correlation kernel is given by c a f,(£, r) = 



TrB[p^B a (t)B b (r)], where the trace is over the environ- 



mental degrees of freedom, B a (t) 



= e lHE 



'■Bt 



iHnt 



, and 



Hb is the Hamiltonian of the environment. Here, p C g is 
the statistical operator of the bath at time 0, which we 
assume to be the one for thermal equilibrium. Hence, 
p°g introduces the temperature into the description of 
the dynamics. This should imply that the system in 
contact with the bath is driven toward a steady state 
that coincides with the thermal equilibrium of the sys- 
tem, p(t — > oo ) oc exp(— /3-ff), with the same temperature 
as the bath [lOj . This is the case if the power spec- 
trum C ab (uj) = J^° dt c a b{t) e~ lult satisfies the detailed- 
balance condition [3| 



Cab 



s /3" 



Cba(oj). 



(3) 



Gaspard and Nagaoka [ll| have shown that the dy- 
namics introduced by the NMME can be obtained not 
only by a numerical integration of Eq. ([T]) but also by 
the solution of a SSE for a state |\l/(£)), 

a 

-t\ 2 Y,S a f dt'c a b{t')e- lflt 'Sb\^{t-t')). (4) 

a,b J ° 

In this non-Markovian SSE (NMSSE) the complex noises 
7a (£) have the properties 



7a(i) = 0, 7a(*)7b(f ) = 0, l*(thb(t') = c ab (t-t') 

(5) 
and one can obtain the dynamics of the open quantum 
system by an average over realizations of the stochastic 
process, indicated by the overline. In particular, the re- 
duced density operator is obtained as p(t) = |\&(£))(\&(t)|. 
However, any attempt to solve the NMSSE (j4]) requires 
a large numerical effort due to the time integral, which 
needs to be evaluated at every time step and for every re- 
alization. This begs the question of whether one can find 
a simpler SSE that reproduces on average the dynamics 
induced by the NMME. In the following we present such 
a SSE. 

We will show that the time-convolutionlcss SSE (TCL- 
SSE) 

id t \*(t)) = (H + \J2la(t)S a -i\ 2 f(t))\*(t)) (6) 

^ a ' 

with 

f(t) = VS a / dt' c ab {t') e~ ikt ' S h e iht ' (7) 

a,b J ° 

reproduces on average the dynamics induced by the 
NMME (HJ) up to fourth order in A. To prove this, 
we write Eq. (J5]) in the interaction picture, |\frj-(£)) = 



e i6t \V(t)) and S a (t) = e i6t S a e- i6t , and expand the 
time-evolution operator up to second order in A, 



!*/(*)} = 



l-*AV / dti-y a (h)S a (h) 

„ Jo 



a . b 



A 2 ^/ dh / dt 2 c ab (t2)S a (t 1 )Sp(ti-t 2 ) 



A 2 V ( dh [ dt 2la 

„ j, JO Jo 



(tl)Sa{tlhb{t 2 )S b (t 2 ) 



x |* 7 (0))+O(A 3 



(8) 



This expansion is inserted into the expression for the re- 
duced density operator pi(t) = \^ j (t)) (^f i (t)\ . By per- 
forming the average, using the properties in Eq. ([S]) and 
the identity c a b(r,t) = c£ a (£, r), and differentiating with 
respect to t we arrive at 



dtPiH) 



„ i Jo 



dT [CabitrfSbWpTMSait) 



c,a(t,T)S a {t)Sb(T)pi(0) 

C: b (t,T)S a (t)pj(0)S b (T) 

c^(t,r)/5 J (0)S 6 (r) ) Sa(t)]+O(A 4 ) 



(9) 



Note that the averages of the terms in A 3 vanish. Further- 
more, replacing pi(0) by pi(r) + 0(X 2 ) docs not change 
the equation up to terms in A 4 . Finally, by returning 
to the Schrodinger picture we arrive at the NMME ([1]) 
up to terms of order A 4 , i.e., higher than the order up 
to which these equations are valid. This is remarkable 
since one might expect a more complex time-non-local 
SSE to be required for reproducing the dynamics of the 
NMME Q. Still, the TCLSSE is local in time, i.e., the 
operator T(t) does not depend on the state of the system 
at previous times and thus can be calculated before the 
numerical integration and be used for each realization 
of the stochastic process. Hence, the numerical cost of 
solving each realization of the TCLSSE is comparable to 
that of a Markovian SSE [4j, 111]. However, the proposed 
TCLSSE requires the generation of colored noise and thus 
the method will only be practicable if an efficient algo- 
rithm for the generation of this noise is available. 

Such an algorithm indeed exists, as we show below. 
We consider only a single bath operator; the generaliza- 
tion to several bath operators is straightforward. Exist- 
ing algorithms for the generation of (real) colored noise 
rely on the numerical solution of a stochastic differential 
equation that has to produce noise with the given target 
correlation function c(t) [12|. However, such an equation 
is a piece of information that is rarely available, since 
even the analytic expression of c(t) may not be known. 
Except for a few simple models, it is more common to 
have access to the power spectrum C{uj). Indeed, C(u) 
is connected with the transitions in the bath. The algo- 
rithm presented in Ref. [13j does not only require C(uj) 



but also the inverse Fourier transform of its square root. 
This quantity is then convoluted with a white noise to 
generate the target real colored noise. We will introduce 
an algorithm that directly uses yJC{uf) as input, thereby 
reducin g th e numerical cost compared to the algorithm 
of Ref. [13[ and that generates a complex colored noise. 
In the following, we present a portable algorithm to 
generate complex colored noise with the properties in Eq. 
([5]), taking a given power spectrum C(w) as input. One 
can easily prove that the noise "f(t) can be generated by 
the relation 



7(*) = 



did 



y/C{w) x(u) e v 



(10) 



where x(ui) is a white-noise process in the frequency do- 
main satisfying 



a;(w) = 0, x(w)x(w') = 0, x*(u)x(uj') = 5(w-w'). (11) 



By substituting the definition (fTU|) into J*(t)j(t') and us- 
ing Eq. (|11[) . we immediately arrive at Eq. ([S])- The other 
relations in Eq. ([5]) arc proven in a similar way. From a 
numerical point of view, the generation of this colored 
noise requires the calculation of the Fourier transform in 
Eq. (fT0|) . A similar algorithm restricted to real noise has 
been proposed in the past 1J, [l5| . 

In order to compare our algorithm with the two from 
Refs. 



13| and 



ompare o 

BE!, 



we choose a test function for 



which wc know c(t) and C(w) analytically, namely c(t) 



-t 2 /2a 



V(2kt 2 ) 1/4 and C(w) = (2ira 2 ) 1/A e" 



/2 



We 



e 

fix a = 1 as our unit of time and choose the inter- 
val t <G [—25,25] for the numerical Fourier transfor- 
mation. To quantify the agreement between the tar- 
get c(t) and the approximations given by the three al- 
gorithms we use the relative error 8 C = J_ dt \c(t) — 

TWwivrL^icwi 2 . 

In Fig. [T] we report the computational cost of calcu- 
lating the correlation function c(t) vs. the number of dis- 
cretization steps of the time interval, TV. Our algorithm is 
moderately faster than the others for any value of N. In 
the inset of Fig.Q] we show the relative error S c calculated 
for 1024 time steps by averaging over an increasing num- 
ber of realizations of the noise. The new method gives 
the best approximation, i.e., the minimum 5 C , for any 
number of runs. While this computational gain is of lit- 
tle importance, the main advantages of the new method 
are that it works directly with the given power spectrum 
and produces a complex noise. 

Finally, we illustrate our results by discussing how the 
TCLSSE is able to reproduce the dynamics induced by 
the NMME and how it describes thermal relaxation for a 
specific system. We consider the coupling of an electronic 
system to the electromagnetic field in a three-dimensional 
cavity. In the dipole approximation, one can derive from 
first principles the power spectrum for this system (we 
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FIG. 1: (Color online) Computation time as a function of the 
number of time steps iV for our method compared to those 
of Refs. [H] and [IlBjJ. Inset: Relative error for N = 1024 
time steps calculated by averaging over an increasing number 
of realizations of the noise. 



set the speed of light to unity) 



a 



cav(w) 



nVe 



[ns(|w|) + 0{—uj)] for \uj\ < lu c 



(12) 

where Ub(uj) = l/(eJ 3u) — 1) is the Bose-Einstein distri- 
bution function, 9{ui) is the Heaviside step function, V 
is the volume of the cavity, and lu c is a cutoff frequency 
determined by the dimensions of the system. This cutoff 
is necessitated by the assumption made in the dipole ap- 
proximation that the electromagnetic field is uniform in 
the region of space of the system. For |w| > uj c , the power 
spectrum is set to vanish |l6|. One can show that the 
detailed-balance condition §3§ is satisfied by this power 
spectrum. Since for this model system the correlation 
function c(t) is not given in an analytical form, we will 
use Eq. ([TU]) to generate the noise. 

In order to quantify the agreement between the noise 
generated by Eq. (|10p and the power spectrum in 
Eq. (|12[) , we have performed a Fourier transform of the 
time-domain signal and compared it with our target. Fig- 
ure H] shows that the agreement is excellent. 

For the electronic system we use a three-site tight- 
binding chain described by the Hamiltonian 



H 



-T {c\c 2 + cjci + c\c3 + C3C2) , 



(13) 



where the operators c] create an electron at site i, and 
assume a single electron to be present. This system is 
coupled to the electromagnetic field inside the cavity by 
the operator 



S = -qJ2^-(Wi\r\W p )cjc p , 



(14) 




FIG. 2: (Color online) Comparison between the target, 
Eq. (|12|l (solid lines) and the Fourier transform into the 
frequency domain of the correlation function obtained from 
Eq. ()10|) by averaging over 90000 realizations of the noise 
(dashed lines). 



where q is the charge of the electron and |Wi) is the state 
localized at site i. For simplicity, we have assumed that 
each mode of the cavity has the same polarization di- 
rection u, parallel to the tight-binding chain. Note that 
the form of this operator should be immaterial for the 
establishment of thermal equilibrium, which is only de- 
termined by the power spectrum. Indeed, we can check 
if the detailed-balance condition is necessary for the sys- 
tem to reach thermal equilibrium. To that end, we use 
the operator in Eq. (|14p within a Markov approximation 
for the correlation function, c(t) oc S(t). We find that a 
steady state is approached that docs not correspond to 
thermal equilibrium [9j. 
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FIG. 3: (Color online) Dynamics of the occupation prob- 
abilities of the eigenstates of the Hamiltonian (|13|) as cal- 
culated from the evolution of the TCLSSE (dashed lines) 
and the NMME (solid lines) with the power spectrum given 
by Eq. (|12|l . The red diamonds represent the thermal- 
equilibrium probabilities calculated from Eq. (I15|) . 

In Fig.[3]we show the probability of occupying the three 
eigenstates of the Hamiltonian as a function of time cal- 
culated using the TCLSSE (dashed lines) and the NMME 
(solid lines), respectively. For the TCLSSE the results 



have been obtained by averaging over 90000 independent 
realizations of the noise. We have used the parameters 
f3 = 1, w c = 1, T = 1, and A = 0.1 and we have used the 
Eulcr algorithm [HI 113 witn time stc P At = °- 005 to 
numerically solve the equations. The dynamics induced 
by the NMME and the TCLSSE are in good agreement: 
The small discrepancies in the numerical solutions are 
due to the finite number of realizations we have used; 
the solution of the TCLSSE still contains some noise, as 
expected. For long times, both formalisms converge to 
the thermal-equilibrium probabilities 



,-Pu 



Pi = 



g-<3ei _|_ e -/3e 2 _|_ e ~' 3e 3 ; 



(15) 



where e, are the eigenvalues of the Hamiltonian (fT3"|) . If 
we were only interested in the long-time limit, we could 
have averaged over all times after some equilibration time 
imin to obtain better statistics, using the ergodic theorem 
to replace the average over many realizations with an 
average over the time evolution of a single realization. 

In conclusion, we have introduced a time-local (time- 
convolutionless) version of a non-Markovian stochastic 
Schrodinger equation, which correctly describes the ap- 
proach to thermal equilibrium as obtained from a non- 
Markovian master equation. This equation can be inte- 
grated with moderate numerical cost, comparable to that 
of a Markovian system. It also shows more advantageous 
scaling with the number of states compared to the master 
equation, which is particularly useful for larger systems. 
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